Ion thruster is one of the electrical thruster types used in spacecraft propulsion and satellites. Electrical thrusters use different acceleration methods for producing the thrust. These methods could be separated in to three classes: electro-thermal, electrostatic, and electromagnetic. Hall thrusters and gridded ion thrusters are in the electrostatic class. In both ion thrusters, ions are electrostatically accelerated out of the plasma. In the gridded ion thrusters, thrust is transferred by the electrostatic force between ions and two grids, but in the Hall thrusters, thrust is created by an external axial electric field established perpendicular to the radial magnetic field.^[1,2] Hall thrusters are well known as the stationary plasma thrusters,^[3] which use noble gases such as xenon, argon, and krypton for propellant. Xenon is usually used because it has high molecular weight and low ionization potential.^[4] The applied potential between the anode and external hollow cathode produces the axial (z-direction) electric field and the outer and inner magnetic polls create the radial (x-direction) magnetic field. Electrons coming from the external hollow cathode are accelerated to the anode by the electric field. However, they will be confined by the radial magnetic field and thus they are trapped in an azimuthal direction (y-direction) around the annular channel. The trapped electrons collide with neutrals and generate ions. The ions with respect to the electrons have larger Larmor radius and therefore they are not affected by the magnetic field. The extracted ions from the system without neutralizing can cause dangerous spacecraft charging. So, it must be neutralized by part of electrons that come from the external cathode. Despite many studies have been done on different aspects of Hall thrusters,^[5–9] there are some problems in physical interpretations of Hall thrusters such as electron anomalous transport across the magnetic field.^[10,11] Among different mechanisms proposed in Refs. [12–14], in order to explain the anomalous electron transport, most physically explanations are relevant to the formation of instabilities in azimuthal, E × B direction. The azimuthal velocity of the ions is almost zero because they are not affected by the magnetic field. However, the electrons are magnetized and they have the azimuthal drift velocity (V_ed = E/B) and the amount is as high as 10⁶ m/s. This difference between the ion and the electron drift velocities could lead to the azimuthal instability with short wavelength and high frequency.^[15,16]

Some researchers have studied the instabilities in the azimuthal direction of the Hall thruster theoretically and experimentally. In 1970, Gary and Sanderson investigated the transfer of magnetized Maxwellian electrons and unmagnetized ions in the plasma. They have shown an instability excited in the azimuthal direction because of coupling between Bernstein and ion acoustic modes.^[17,18] In 2004, Litvak et al.^[19] experimentally identified and characterized the high-frequency plasma azimuthal waves in Hall thrusters. In 2012, McDonald^[20] experimentally showed that the azimuthal drift instability could increase the electron mobility in the axial direction of the thruster. Lafleur et al.^[21] have used one-dimensional (1D) azimuthal particle in cell (PIC) simulation to study the effect of the plasma density on the electron transport perpendicular to the magnetic field in Hall thrusters. They found that the electron transport can be explained by the classical electron–neutral collision theory in low-density plasma. However, in high-density plasma, a strong instability can be observed which causes considerable increase in electron transport mobility, even without considering the electron–neutral collisions with respect to the classical theory.

Although simulation in two or three dimensions is more accurate comparing to that in one dimension but it requires long time. Therefore, it is useful to study the simulation in one dimension in the azimuthal direction to investigate the azimuthal instability.

In this work, the azimuthal E × B drift instability is studied in the Hall thruster using one dimensional PIC and Monte Carlo collision (MCC) simulation method. In Section 2, the theory of E × B drift instability is studied by using the dispersion relation in the Hall thruster and it is used to investigate the effect of the electric and magnetic fields on the number of unstable modes. PIC simulation and detailed conditions are briefly explained in Section 3. In Section 4, the electric field in the azimuthal direction and the azimuthal drift instability are characterized. Furthermore, the effect of the electric and magnetic fields on the azimuthal instability frequency and wavelength is studied. The profiles of electron and ion velocity distribution functions and the phase space diagram are discussed in Section 5. Finally, in Section 6, the results are summarized.

Many oscillation modes in the range of a few up to several tens of MHz^[22] have been identified in Hall thrusters. These modes propagate in radial, azimuthal, longitudinal, and total (sum of the above three components) directions. One of the well-known oscillations in azimuthal direction is the electron drift instability which arises from coupling between Bernstein and ion acoustic modes.

In order to understand the origin of the E × B drift instability, we need to study the dispersion relation. In this paper, the dispersion relation is studied using the Gordeev function because in this method, the solution of the dispersion relation converges very quickly.^[23] Considering constant electric and magnetic fields, where the axial electric field is perpendicular to the radial magnetic field, the dispersion relation is obtained by the electrostatic growth of the wave in uniform plasma in the Hall thruster. The dispersion equation in the Hall thruster is^[16,23] 1a 1b where g(Ω, X, Y) is the Gordeev function, λ_De is the electron Debye length, k_x, k_y, and k_z are the x, y, and z components of the wave vector k, respectively, ω_pi, ω_ce, and ω are the ion plasma, electron cyclotron, and plasma frequencies, respectively. In addition, ρ_e = v_th/ω_ce is the electron Larmor radius, is the electron thermal velocity, k_B is the Boltzmann constant, T_e and m_e are the temperature and mass of electrons, respectively. Moreover, Z(z) and I_m(x) are the plasma dispersion function and the modified Bessel function of the first kind, respectively. With normalizing the lengths to λ_De, the frequencies to ω_pi, and the velocities to C_s, where C_s is the ion sound speed, the normalized dispersion relation is given by 2 where is the ratio of the xenon mass (m_i = 2.18× 10⁻²⁵ kg) to the electron mass. Under one-dimensional approximation k = (0, k_y, 0), the linear dispersion relation shows that the unstable modes only occur in separated lines corresponding to the electron cyclotron frequency harmonics^[24] 3 where n is an integer. By solving Eq. (2) using the Gordeev function and similar method presented in Ref. [14], eight unstable modes corresponding to the cyclotron frequency harmonics are found, as shown in Fig. 1. In addition, k_y in Eq. (3) becomes 4 Since λ_De = 7.43 × 10⁻⁵ m, the wavelength corresponding to the first harmonic of the electron cyclotron frequency is λ₁ = 1.76 × 10⁻⁵ m.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. The effect of kx on (a) angular frequency (real part of ω) and (b) growth rate (imaginary part of ω) at kz = 0. The used parameters are presented in Table 1.

Table 1.

Table 1.

Table 1. The wavelength of the first harmonic of the electron cyclotron frequency for different Bx. . B/T 0.01 0.012 0.014 0.016 0.018 0.02 0.022 0.025 λ1/mm 7.14 4.96 3.64 2.79 2.2 1.76 1.47 1.14

Table 1. The wavelength of the first harmonic of the electron cyclotron frequency for different Bx. .

It is obvious from Fig. 1 that for k_yλ_De ≈ 0, there are delta-function like peaks. However, since k_yλ_De = 0, it is concluded that g(Ω, 0, 0) limits to infinity. Therefore, the solution of the dispersion relation is infinity and we have the resonance condition. The real part of ω (angular frequency) is shown in Fig. 1(a). Figure 1(b) is also related to the imaginary part of ω known as the growth rate of the instability. One can see from Figs. 1(a) and 1(b) that by increasing , the peak height decreases which leads to disappearing of both growth rate and angular frequency curves.

Furthermore, the effect of the electric and magnetic fields on the number of unstable modes is studied. In Fig. 2, the effect of the magnetic field on the number of unstable modes is depicted. It is found that by increasing the radial magnetic field (with constant electric field E_z = 2000 V/m), the number of unstable modes decreases approximately proportional to 1/B². It is concluded from this figure that in order to have a stable mode (for fixed electric field), a larger magnetic field can be useful. Figure 3 shows the unstable modes as a function of electric field (with constant magnetic field B_x = 0.02 T). This profile is approximately linear. As it is seen, the system is unstable for a large axial electric field. The physical interpretation is that when the axial electric field is high, more energy transfers to the fluctuations of the system and consequently the unstable modes increase.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. The effect of magnetic field on the unstable modes with constant electric field E0 = 2000 V/m.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. The number of unstable modes as a function of the electric field with constant magnetic field B0 = 0.02 T.

In Table 1, the wavelength corresponding to the first harmonic of the electron cyclotron frequency for different radial magnetic fields is shown. In Table 2, the wavelength corresponding to the first harmonic of the electron cyclotron frequency for different axial electric fields is shown.

Table 2.

Table 2.

Table 2. The wavelength of the first harmonic of the electron cyclotron frequency for different Ez. . Ez/V·m−1 5000 10000 12000 14000 16000 18000 20000 22000 25000 30000 λ1/mm 0.44 0.88 1.07 1.25 1.43 1.61 1.76 1.14 2.23 2.67

Table 2. The wavelength of the first harmonic of the electron cyclotron frequency for different Ez. .

The results could be physically interpreted as follows. By increasing the magnetic field, more electrons are confined and thus the number of unstable modes and consequently the wavelength decreases, while by increasing the electric field, the electron drift velocity increases to distort the electron confinement, which leads to more unstable modes with longer wavelength. Although a small electric field has less unstable modes and smaller amplitude of unstable modes but it is not good for Hall thruster because of the low ion beam velocity.

In Fig. 4, the effect of different ions on the unstable modes is investigated. As it is expected, the ion mass does not affect the number of unstable mode because in Eq. (3) there is no ion mass term. However, it can affect the amplitude of the unstable modes. Even Xe leads to the high amplitude instability, but it has high molecular weight and low ionization potential, then it is suitable as a propellant.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Effect of ion mass on the (a) angular frequency and (b) growth rate at kz = 0.

In order to obtain the velocity distribution of electrons, it is useful to study the plasma system using the kinetic theory. In other words, in the kinetic approximation, the plasma system is characterized by the electron distribution function. The analytical solution of the kinetic equation is complicated and therefore numerical solutions such as Vlasov^[25] and PIC^[26–28] simulation methods should be used. In PIC simulation, the system is described by the electromagnetic fields, distribution functions, and Newton’s equations.

In this paper, the numerical method is organized by 1D azimuthal-PIC-MCC simulation algorithm.^[29] Radial magnetic field B₀ and axial electric field E₀ are uniform and fixed during the simulation. At the beginning of the simulation, the charged particles (ions and electrons) are separated uniformly in the axial (acceleration length) and azimuthal domain of the simulation. Both initial velocity distribution functions of electrons and ions are Maxwellian. All physical parameters used in this work are summarized in Table 3, in which L_y, L_ch, and L_{z, acc} are the azimuthal, axial, and axial acceleration lengths of the thruster channel, n_n and n₀ are the plasma and gas densities, T_n is the temperature of neutral gas, and dt is the time step in the simulation.

Table 3.

Table 3.

Table 3. Physical parameters used in the PIC simulation. . Parameter Value Ly/m 0.0267 Lz/m 0.025 Lz, acc/m 0.010 ne/m−3 1 × 1017 n0/m−3 1 × 1018 B0/G 200 E0/(V/m) 10000 Ti/K 2320 Te/K 116000 Tn/K 700 dt/s 2 × 10−12 Nc 500 Npcc 500

Table 3. Physical parameters used in the PIC simulation. .

The self-consistent electric field is solved along the E × B direction using Thomas tridiagonal algorithm with periodic boundary conditions.^[30] The reference potential at the boundaries is zero. Furthermore, in this simulation, the axial boundary conditions are considered. The electrons and ions exited from the axial region are replaced by new electrons with temperature 10 eV and ions with temperature 0.2 eV (see Fig. 5). In addition, the electrons are magnetized and the electron–neutral collision is the elastic one. The neutral atoms are assumed to be once ionized, and during the ionization processes, only the energy loss and scattering are assumed. Therefore, new electrons and ions are not included in the simulation. Moreover, it is supposed that the ions are unmagnetized and the ion–neutral collision frequency is neglected.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. Schematic of simulation domain used in the 1D PIC simulation.

For checking the code, the average kinetic energy of electrons and ions is plotted in Fig. 6. It is clear from this figure that the results have a good behavior since the electron kinetic energy reaches to the steady state at t > 2 μs. This indicates that the code can be useful for explanting the physical phenomena. Since heavy ions are accelerated in the axial direction, their kinetic energy is higher than that of electrons. The electron and ion density fluctuations at t = 12 μs is shown in Fig 7. Calculating the average electron and ion fluctuations over last 10 μs shows that the electron fluctuation is approximately 20% while the ion density fluctuation is approximately 40%. In this simulation, N_pcc = 500 cells and N_c = 500 particles per each cell, and thus Δy/λ_De = 0.75.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. Mean energy of electron and ion with number of particles per cell Npcc = 500, number of cells Nc = 500, and azimuthal length Ly = 0.026 m. The parameters are presented in Table 1.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. Electron and ion density fluctuations at t = 12 μs. The used parameters are presented in Table 1.

3.1. Electric field in azimuthal direction

As it is known, the azimuthal instabilities could lead to non-uniformity and polarization of the plasma in the azimuthal direction. These instabilities are excited by the formation of an azimuthal electric field. However, the electron confinement can be destroyed by the azimuthal electric field. The oscillating electric field in the azimuthal direction as a function of position (y) at t = 3 μs is depicted in Fig. 8(a) and the time evolution of the electric field at constant position (y = 1.78) cm is shown in Fig. 8(b). However, this azimuthal electric field can be correlated with the plasma density (Fig. 6) in the Hall thruster and generate a pure axial electron current.^[31] Therefore, the efficiency of the Hall thruster decreases.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. (a) Oscillation of the azimuthal electric field at t = 3 μs and (b) time evolution of the azimuthal electric field at y = 0.018 m. (c) Time–space evolution of the azimuthal electric field up to t = 8 μs and Ly = 0.026 m The used parameters are presented in Table 1.

The wavelength, frequency, and velocity of the E × B drift instability can be characterized by the contour plot (time and space evolution) of the azimuthal electric field which is given in Fig. 8(c). It is clear that after a short time, an obvious monochromatic wave can be seen. The wavelength of this instability can be obtained by counting the numbers of modes per length in y-direction that is 1.67 mm and the instability frequency can be obtained by counting the number of modes per time which is 4.5 MHz. Thus, the propagation velocity of the azimuthal instability is 7650 m/s. The obtained results are in good agreement with the reported one in Refs. [23,32,33].

Since particles motion depends on the electric and magnetic fields, it is useful to investigate the effect of the magnetic and electric fields on the propagation velocity of the instability. In Table 5, the effect of the magnetic field on the wavelength, frequency, and velocity of the instability has been illustrated. It is clear that by increasing the magnetic field, the wavelength of the instability decreases while the instability growth rate increases. Finally, by increasing the magnetic field, the instability velocity decreases. Since the Hall current^[1] (I_H ∝ V_d/B, where V_d is the discharge voltage) decreases by increasing the magnetic field, it is expected that the instability velocity decreases. Furthermore, by increasing the magnetic field, the electron drift velocity V_ed decreases and consequently the difference between the electron and ion velocities decreases, which leads to the decrease of the instability velocity. Comparing Table 1 and Table 4, we obtain a good agreement between the theory and PIC simulation for the wavelength of the instability.

Table 4.

Table 4.

Table 4. The characteristics of instability with different magnetic fields. . B/T 0.01 0.012 0.014 0.016 0.018 0.02 0.022 0.025 λpic/mm 3.8 3.4 3.0 2.7 2.0 1.7 1.5 1.3 ωpic/MHz 3.0 3.2 3.5 3.8 4.2 4.5 4.8 4.9 vφ/m·s−1 11400 1088 10500 10260 8400 7650 7200 6370

Table 4. The characteristics of instability with different magnetic fields. .

Table 5.

Table 5.

Table 5. The characteristics of instability with different electric fields. . Ez/V·m−1 5000 10000 12000 14000 16000 18000 20000 22000 25000 30000 λpic/mm 0.8 0.9 1.0 1.2 1.4 1.6 1.7 1.9 2.1 2.4 ωpic/MHz 5.6 5.5 5.0 4.75 4.62 4.5 4.5 4.25 4.0 3.75 vφ/m·s−1 4480 4950 5000 5700 6468 7200 7650 8.075 8400 9000

Table 5. The characteristics of instability with different electric fields. .

It is shown in Table 5 that by increasing the electric field, the wavelength of the instability increases while the instability frequency decreases. Furthermore, by increasing the electric field, the velocity of the instability increases. The physical interpretation of these behaviors can be given as follows. By increasing the electric field, the electron drift velocity V_ed increases and consequently the difference between the electron and ion velocities increases, which lead to the increase of the instability velocity. Moreover, increasing the electric field needs to increase the discharge voltage that leads to the increase of the Hall current, consequently the velocity of the instability increases. By comparing between quantities presented in Table 2, one can see that the results obtained by the theory are in good agreement with the ones obtained by PIC simulation for the wavelength of the instability.

Finally, one can conclude that the dependence of the wavelength, frequency, and velocity of the instability on the magnetic field is more than their dependence on the electric field.

3.2. Velocity distribution functions and phase space

The velocity of particles plays a key role in the interaction between travelling waves and charged particles. Therefore, it is necessary to distinguish between the velocity distribution functions of charged particles. The kinetic equations give us a suitable explanation of particle–wave interaction. Using these equations, the electron and ion densities and velocity distribution functions f_i,e(r v t) can be obtained. It is important to note that the distribution functions of electrons and ions play a key role in predicting and controlling the plasma parameters.

The correlation between the oscillating electric field and distribution function has a strong effect on the distribution function shape at different time. In Fig. 9(a), the profile of the electrons velocity distribution function (EVDF) at five different time is plotted. One can see from this figure that as the time passes, the width of the electron distribution function increases, which means that the electrons temperature increases. In other words, the number of fast electrons increases. The relationship between the velocity and temperature of the electrons can be obtained by the partition function 5 In Fig. 9(b), the electron temperature as a function of y is depicted at different time. By using Eq. (5) and calculating the average electrons temperature in all azimuthal lengths, it is found that at t = 0.5 μs the electron average temperature is 30 eV. Moreover, at t = 3 μs the electron average temperature is 64.8 eV, being more than twice. In addition, the electron average temperature at t = 11 μs is 65 eV, which is close to that at t = 3 μs. This is good because the width of the electron distribution function at t = 3 μs and t = 11 μs is approximately the same.

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. (a) Electron velocity distribution function and (b) electrons temperature in azimuthal length at different time.

In Fig. 10(a), the distribution function of ions is plotted for different time. One could see that with the passage of time, the ion velocity distribution function deviates from the Maxwellian distribution function. Somewhere, there are approximately two peaks indicating the ions drift in the azimuthal direction. This confirms that there is instability in the system. We also superpose in Fig. 10(b) a best fit using a double-Gaussian model at t = 6 μs according to 6 The empirical fit of the simulated distribution was first used in Refs. [34–36].

Fig. 10.

	Figure Option ViewDownloadNew Window
	Fig. 10. (a) Ion velocity distribution function, (b) fitted curve at t = 6 μs, (c) ion temperature in azimuthal length at different time.

The ion temperature in azimuthal direction for three different time is plotted in Fig. 10(c). By getting average over all azimuthal lengths and comparing the results, one can conclude that by increasing the time, the width of the distribution function increases. This confirms that the average temperature increases. Another important point that can be obtained from this figure is that the average ion temperature at t = 0.5 μs is 0.225 eV. At later time, i.e., at t = 3 μs and t = 11 μs, the average ion temperature is 3.23 eV and 3.05 eV, respectively, which is approximately constant. These results illustrate that the average ion temperature at t = 3 μs and t = 11 μs is more than ten times of that at t = 0.5 μs.

It is important to note that the main reason for the deviation of the distribution function from the initial state is the nonlinear process. However these lead to the saturation of the instability. Usually, the saturation process occurs after increasing the instability growth rate up to several times. The saturation happens because of ion–wave trapping. The study of ions and electrons phase space at two different azimuthal and axial directions shows that in the axial direction, there is no evidence of trapping for both electrons and ions. Furthermore, the phase space diagram in the azimuthal direction for electrons shows that there is no electron trapping, but for ions, phase space vortexes are created as illustrated in Fig. 11. Similar ion trapping was reported in many other simulations where the instability was investigated.^[37,38]

Fig. 11.

	Figure Option ViewDownloadNew Window
	Fig. 11. Ion phase space in azimuthal direction of Hall thruster at t = 3 μs.

The ion mass is one of the most important parameters in Hall thruster operation. Considering different cross section of excitation and ionization energy for each gas (Table 6), we investigate the effect of the ion mass on the instability. The used parameters in Table 3 expect L_y = 1.34 cm.

Table 6.

Table 6.

Table 6. Ionization and excitation energies for noble gases. . Gas Mi/kg Ionization energy/J Excitation energy/J He 6.64 × 10−27 3.938 × 10−18 3.1723347 × 10−18 Ne 3.35 × 10−26 3.455 × 10−18 2.6468166 × 10−18 Ar 6.63 × 10−26 2.524 × 10−18 1.8505286 × 10−18 Kr 1.39 × 10−25 2.243 × 10−18 1.5877695 × 10−18 Xe 2.18 × 10−25 1.943 × 10−18 1.3330212 × 10−18

Table 6. Ionization and excitation energies for noble gases. .

By using the contour plot of the azimuthal electric field for different ion masses (Figs. 12(a)–12(e)), the wavelength, frequency, and phase velocity of the instability are shown in Table 7. It is clear that the ion mass cannot affect the instability wavelength, which is good agreement with the theory (see Fig. 4(a)). But the instability frequency and phase velocity extremely depends on the ion mass. The instability characteristics for Ar, Kr, and Xe are also in good agreement with the E × B electron drift instability. While for light mass ions such as He and Ne, the frequency is out of range for drift instability frequency (1–10 MHz).

Fig. 12.

	Figure Option ViewDownloadNew Window
	Fig. 12. Time–space evolution of the azimuthal electric field up to t = 12 μs for different ions: (a) He, (b) Ne, (c) Ar, (d) Kr, and (e) Xe.

Table 7.

Table 7.

Table 7. Characteristics of instability for noble gases. . Gas λpic/mm ω/MHz vφ/m·s−1 He 1.67 25 41750 Ne 1.67 13.5 22545 Ar 1.67 8 13360 Kr 1.67 5 8350 Xe 1.67 4.9 8183

Table 7. Characteristics of instability for noble gases. .

The dispersion relation was investigated using Gordeev function to understand the drift instability theory in Hall thruster. The results indicated that the unstable modes occur only in the separated lines corresponding to the electron cyclotron frequency harmonics. Furthermore, the effect of the electric field, magnetic fields, and ion mass on the number of unstable modes was studied. It was found that the increasing axial electric field leads to the increase of the number of unstable modes. While by increasing the radial magnetic field, the number of unstable modes decreases and ion mass does not affect the number of unstable modes. Furthermore, by using 1D Azimuthal PIC-MCC simulation, the drift instability was developed in the Hall thruster. The results confirmed that there is an instability with short wavelength in the range of millimeter and high frequency in the range of megahertz. Moreover, the results showed that the wavelength, frequency, and velocity of the instability depend more on the magnetic field than the electric one. Also, the result showed that for light mass ions, the frequency and phase velocity are very high that lead to higher electron mobility in the axial direction.

In addition, the effect of the drift instability on the electron and ion distribution functions was investigated using 1D PIC simulation. The results indicated that in the presence of the azimuthal instability, the electron and ion distribution functions deviate from the initial state (Maxwellian) in the azimuthal direction. Finally, ion trapping in the azimuthal direction confirmed the instability saturation process.